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'"^ ' Abstract. We present dynamical distribution functions for a homogeneous sample of oxygen-rich, evolved, 
^^O , intermediate-mass stars in the inner galactic plane. We use an axisymmetric, two-component Stackel potential that 
^ ' satisfies recent constraints on the galactic potential, amongst others a slightly declining local rotation curve. We show 
^ , that this potential is adequate to model stellar-kinematic samples with radial extent ranging from ~ 100 pc to ~ 5 
kpc in the Galaxy. 

1 The stable two-integral model that gives the best fit to the first three projected moments provides a very good global 
' representation of the data but fails to reproduce the central line-of-sight dispersion, the central apparent scaleheight 
and the almost-cylindrical rotation at intermediate longitudes (5°< \£\ < 15°). All these features, indicative of the 
galactic Bar, are fitted well by a three-integral model. We discuss various properties of the two- and three-integral 
, distribution functions and the implications for galactic structure. A somewhat thicker disk component is needed to 
' explain the observed distribution of older AGB stars in the plane; this component at the same time fits the kinematics 
of AGB stars at higher latitudes better than the thinner disk. We find that the Disk and the Bulge, as traced by AGB 
I stars, are very similar dynamically and could well be one and the same component. There is a dynamically distinct 
2 component in the inner 100 pc of the Bulge, however. 

t2 

> ■ 

■ 1—1,1. Introduction 



X 



^ ■ Methods to analyse observational data almost always fall into one of two classes, "direct" and "indirect" methods. 
5^ , The former seek to derive (deproject) the desired quantities in a direct manner from the observed quantities, the 
latter to predict the observables from a purely theoretical model and then accept or reject the model by comparing 
the predictions to the observations. 

In the field of modelling galactic stellar dynamics, two types of indirect methods prevail. One is to construct models 
via N-body simulation (eg. Fux 1997), the other via the (semi-direct) Schwarzschild method (Schwarzschild 1979, eg. 
Zhao 1996). In this paper we will use an indirect, Schwarzschild-type method to model the stellar dynamics of the inner 
Milky Way Galaxy. We test assumed dynamical distribution functions for their ability to reproduce the distribution of 
our sample of evolved, intermediate-mass stars (Sevenster et al. 1997a, b, S97A,S97B). This sample is representative of a 
large fraction of the stellar content of the Galaxy, but does not sample the old, spherical Bulge or the Halo. Therefore, 
it is justified to consider its dynamical distribution in a global, fixed potential, unlike most Schwarzschild models. 
Rather than trying to build self-consistent models from an unsuitable sample, we constrain the model gravitational 
potential with a variety of other recent observations. The goal is to find the dynamical characteristics of the inner 
Galaxy and whether there are clearly distinct dynamical components. 

In Sect. 2 and Sect. 3 we describe the method and its detailed implementation, in particular the choice of the 
potential. In Sect. 4 we discuss the resulting two-integral model and its errors and stability. We present a three- 
integral model in Sect. 5 and a two-integral model for a galactic-centre sample in Sect. 6. We interpret the results in 
Sect. 7 and we end with conclusions in Sect. 8. 
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Fig. 1. Cuts in longitude and latitude through the OH/IR data smoothed with round kernels (initial kernel 1°, solid), 
smoothed with elongated kernels (initial kernel 1° in £ and 0?5 in 6 , dotted) and the COBE-DIRBE surface-density 
map (dashed). All densities are normalized separately to have a peak density of 1 and interpolated onto a grid of (1°, 
0°67 ). 



2. Method 

The distribution function of a stellar system is a function of at most three isolating integrals of motion Ji , according 
to Jeans' theorems. It gives the density of stars in the full six-dimensional phase-space (x,V). Integrating over all 
velocities, we get the true n*^-order moments M^"^ of the distribution function according to : 

n r ^ 

M(")=p((nvt))= /(n^k)/(I)d3V n = 0,l,... VkG (Vx,Fy,V;) 1 

k=l k=l 

leaving out the dependencies on (x). There are one zcroth -order moment (M^*^)), three first-order moments (M'^'^^), 
nine second-order moments (M^^^) (six of which are independent) and so on. The method we use to model the 
distribution function of a galactic stellar sample was developed by Dejonghe (1989). For details we refer to that 
article; here we discuss the method only briefly. In a given gravitational potential, a distribution function is built from 
a library of orbital components that are (analytic) functions of the integrals of motion in that potential. By minimizing 
the quadratic differences D between moments of the model distribution (Mm) and the observed distribution (Mq), 
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it determines, sequentially, the best combination of components and the corresponding coefficients. The moments can 
have different weights Wi in the determination of D according to the importance they should have in the fit. Because 
of its quadratic-programming character we will use "QP" to refer to the modelling program. There is no true 
connected to the fit, because there is no optimization of free parameters in the strict sense: the parameter D can be 
used only to compare the goodness of fit between models with the same potential and input data. To compare different 
potentials or data sets, the ratio of the initial to the converged value of D might be used. 



3. Implementation 

In this paper, we use QP with an axisymmetric potential. The Galaxy's density distribution is not axisymmetric, but 
the probably small eccentricity of the potential and the not too-strongly barlike inner stellar kinematics (see Sevenster 
1999) indicate that the non axisymmetric part of the potential is negligible in a first approach. The infiuence of the 
third integral is not negligible, in any case not for the galactic Disk (eg. Oort 1965). After starting our investigations 
with two-integral (21) models, that are easier to interpret, we construct a three-integral (31) axisymmetric model to 
try and overcome the limitations of the 21 model. 

In the models presented in this paper, we include the first three projected moments of distribution functions in 
the fit. The moments used for the comparison between model distribution function and observations are hence S (the 
surface-number density), I](Vios) and E(V|^g) . The weights Wi in Eq. 2 are all equal to 1 in the models presented in 
this paper. In the figures we will mostly show the more commonly-used derived moments S, (Vlos) and aios- 
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Fig. 2. The longitude-velocity diagram for the OH/IR stars with the rotation curves for the BD2 (long-dashed), the 
HI2 (short-dashed) and the AX2 (solid) potentials (see Sect. 3.2, Table 2). The thin dashed hnes connect the maximum 
and minimum observed velocities, respectively, in longitude bins, stellar velocities. The BD2 potential yields too high 

velocities at intermediate longitudes; since stars have dispersions the rotation curve should fall in general somewhat 
below the maximum observed stellar line-of-sight velocities, as is the case for the AX2 and HI2 potentials. 

Table 1. Average dispersions in coordinate separations. 
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3.1. Data 

The data were acquired specifically to constrain optimally dynamical models of the galactic Plane (S97A, S97B). 
The sample consists of positions on the sky (accuracy ~0'.'5) and line-of-sight velocities (accuracy ~1 kms~^ ) with 
respect to the local standard of rest (LSR) of OH/IR stars; oxygen-rich, asymptotic-giant-branch (AGE) stars in the 
thermally-pulsing phase. These stars form a partly relaxed population (0.5-7.5 Gyr, Sevenster 1999) and trace the 
dominant mass distribution (Frogel 1988). The region covered in galactic coordinates is —45° < £ < 10° and |6| < 3°. 
In total 507 objects were found, forming the AOSP (Australia telescope Ohir Survey of the Plane) sample used in this 
paper. The QP program will correct the velocities for the motion of the LSR, assuming the LSR is on a circular orbit 
at Rq = 8 kpc in the model potential. 

In Table 1 we give the dispersions, using all stars in the sample, in the distribution of separations in all three 
coordinates for different numbers of nearest neighbours TVnn (on the sky). The average velocity difference between 
stars does not change with number of nearest neighbours. Since the velocity profile sampled by the stars has to change 
with position on the sky, this means that the velocities of neighbouring stars are completely independent. Therefore, 
we use adaptive-kernel smoothing to grid the data on the sky (Merritt & Tremblay 1994), but treat the velocity 
coordinate separately (cf. their equation 40). 



4 Sevenster et al.: Distribution functions for evolved stars 



Table 2. Potential models. Rq = 8 kpc. This table is discussed in Sect. 3.2. 
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Fig. 3. Circular-velocity curves for various model potentials and as derived from observations. In both panels the 
thin lines indicate observed rotation curves, derived from HI data (solid, Burton & Gordon 1978) and from the stellar 
sm+acc brightness (dot dashed, = 184 km s^^ ^kp^; Allen et al. 1983). The thick curves arc the fitted potentials, 
a Two-component Stackel KK-model potentials fitted to HI (HI2, dashed), to the stellar surface brightness in the 
inner regions (AX2, dot-dashed) and to constant circular velocity outside 5 kpc (BD2, Batsleer & Dejonghe 1994). b 
Three component Stackel KK-model potentials fitted to HI (HI3, dashed) and to the stellar surface brightness (AX3, 
dot-dashed). 



First, the data arc smoothed with initial gaussian kernels of l°xl°x30 kms""'^ . (These initial-kernel sizes were 
optimized to retain the scales of the large-scale distribution without showing individual stars, cf. Table 1.) Then, for 
each star, the spatial kernel is adapted according to the surface density and the mean velocity and velocity dispersion 
are determined from the velocity profile thus created at its position on the sky (so, the ratio of the spatial-kernel sizes 
was kept constant). Finally, the surface density, mean velocity and velocity dispersion are calculated on a regular grid, 
still using gaussian distributions in all three dimensions, with these final parameters. 

In Fig. la-d we compare the resulting surface density with that of the COBE-DIRBE observations (Dwck et 
al. 1995). In the same figure, we show the surface density resulting from smoothing the data with elongated kernels, 
to reflect the possible difference between the vertical and the radial density scale. The COBE- and the AOSP surface 
densities clearly trace a similar population (the evolved late-type stars) and The round kernels provide slightly better 
agreement with the COBE data and are also favoured by the results give in Table 1. We thus use the round-kernel 
surface density for our standard model, but also give results for the elongated-kernel surface density. 
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Table 3. Observed values for the potential parameters (units as in Table 2). 



Quantity 


Value (rcf) 


Value (ref) 


Value (ref) 


Value (ref) References 


A 


14.4 ± 1.2 (1) 


11.3 ± 1.3 (2) 


19 ± 6 (3) 


14.82 ± 0.84 (9) 1 Kerr & Lynden-Bell 1986 


B 


-12.0 ±2.8 (1) 


-13.9 ±0.9 (2) 


-13 ±5 (3) 


-12.37 ±0.64 (9) 2 Hanson 1987 


Vlsr 


184 (4) 


200 ± 10 (10) 231 ± 21 (9) 


3 Evans & Irwin 1995 


2ARq 


228 (5) 


248 (6) 


257 (7) 


252 (9) 4 Rohlfs et al. 1986 


dV/dR 


-3.7 (4) 


-2.4 (9) 




5 Caldwell & Coulson 1989 


Po 


0.1 (8) 


0.076 ± 0.015(11) 




6 Schechter et al. 1989 


«o 


36 (1) 






7 Pont et al. 1994 




46 ±9(12) 






8 Kuijken & Gilmore 1989b 


Af50kpc 


4.9±1.110"(13) 






9 Feast & Whitclock 1997 

10 Merrificld 1992 

11 Crczc ct al. 1998 

12 Kuijken & Gilmore 1989a 

13 Kochanek 1996 



3.2. Potential 

To model the galactic potential, we use so-called Stackel (S) potentials (see de Zeeuw 1985), because for those three 
integrals of motion are known analytically. The specific form of the S-potentials used in this work is that of a multiple, 
axisymmetric Kuzmin-Kutuzov (KK) potential (see Dcjonghe & de Zeeuw 1988). The Galaxy is thus simulated by a 
small number of separate, axisymmetric components with different flattenings, but the same focal lengths to keep the 
over-all potential in Stackel form, ie. separable in ellipsoidal coordinates. Batslecr & Dejonghe (1994) give values for 
the parameters that optimize the fit of a double S-KK potential to the large-scale rotation curve of the Galaxy, with 
the two components representing a (dark) halo and a disk (BD2, Table 2). 

In Fig. 2, we show the longitude-velocity diagram of the AOSP sample, together with the rotation curves for the 
two-component potentials from Table 2. The BD2 potential has little mass in the central regions of the Galaxy; the 
rotation curve is shallower than the i?*^'^ curve determined for the inner Galaxy (Allen et al. 1983; see Fig. 3 and 
the AX2 potential in Fig. 2). The circular velocity at intermediate longitudes (10° to 20°) is too high; the extreme 
stellar velocities should be somewhat larger than the circular velocity, give or take statistical fluctuations, because the 
line-of-sight dispersion is larger than the asymmetric drift (^ o"^/120kms^^ in the Disk). In test runs we found that 
these short-comings inhibit the construction of realistic models for the AOSP sample, that is dominated by the Bulge 
potential. We therefore tried to find an S-KK potential that has more mass in the central regions, yields a realistic 
rotation curve and gives acceptable values for important parameters such as the local circular velocity and the Oort 
constants. 

In Table 2, we list the values for these parameters for three double and two triple S-KK potentials. These potentials 
were constructed to fit various rotation curves for the Galaxy (Fig. 3). The columns in Table 2 give the total mass of 
the Galaxy Mtot, Oort's constants A and B, the local circular velocity Vlsr? the common combination 2ARq, the first 
radial derivative of the circular velocity dV/dR {Rq dV/dR = I^lsr ~ 2A_R0), the local density Pq , surface density 
Eq and local epicyclic frequency Kq, the flattening of the halo ghaioj the flattening of the second (third) component 
gd,b; the fraction of the mass in the second (third) component /d,b and the square of the "focal length" A^. For all 
details on these S-KK potentials and the parameters see Batsleer & Dejonghe (1994). 

In Table 3 we list observed values for some of the quantities in Table 2, determined by various authors. The total 
mass of the Galaxy, the mass fractions of the Disk and Bulge and the flattenings of the components are not well 
established observationally and treated as mere parameters for the potentials rather than physical quantities. Even 
for the double KK potentials there is considerable freedom to create rotation curves of all sorts and at the same time 
obtain very realistic values for the important parameters. HI2 and HIS are based on the assumption that the HI gas 
follows purely circular orbits, which is probably not the case in the inner Galaxy. They are therefore not likely to be 
realistic, but it is interesting that the Hl-rotation curve as well as local parameters can be reproduced with so simple 
a potential. 

An obvious limitation is that A has to be the same for all components, in order to keep the total potential in 
Stackel form. This means that all components simultaneously become more compact when increasing the contribution 
of the inner regions of the Galaxy. Therefore, all five potentials listed in Table 2 give quite acceptable values for all 
observational constants except pq. Accordingly, the vertical forces at larger radii (R ^ Rq) are not in agreement with 
observations (eg. Kuijken & Gilmore 1989a). Related results should be viewed with care; we will give most attention to 
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the inner Galaxy, R ^5 kpc, where the potential is realistic and also the observations sample the distribution optimally 
(Fig. 5) 

We will use the potential AX2, for which we deem the rotation curve most realistic. Its kinematic parameters, 
especially kq, Vlsr- 2Ai?Q, dV/dR and B, arc best in agreement with observations. AX2 docs not have a constant 
outer rotation curve, consistent with recent claims (Rohlfs et al. 1986; Binney & Dehnen 1997; Feast & Whitelock 1997; 
Honma & Sofue 1997). The self-consistent density for the AX2 potential is positive everywhere. Its central scalelength 
(200 pc) and scaleheight (120 pc) are of the order of those found for the density distribution of the AOSP sample 
(Sevenster 1999). 

3.3. Two-integral orbital components 

We use two families of 21 orbital components - fuU distribution functions in themselves - to build the total distribution 
function; the first has infinite extent ("bulge-like"), the second is limited in the vertical direction ("disky"). Their 
functional forms are: 

Fl{a,l3,is) = E" {ELl/2f for isL^>Oand Fl = otherwise , 3 

F2{a,l3,-i,za,is) = S" {2SLlf [{E - S)l{So- S)Y' for > and F2 = otherwise, 4 

with E the total energy and Lz the angular momentum, the two classical integrals in axisymmetric systems. S is the 
energy of an orbit that reaches just to zq out of the plane and So the energy of a circular orbit in the plane, both for a 
given L^. We will discuss the properties of Fl and F2 briefly and refer to Batslcer & Dcjonghe (1995) for a thorough 
treatment of components of these types. The Fl and F2 components are all even in L^; to get rotation, the parameter 
is is introduced. If = — 1 only the co-rotating half of phase space is populated, for ig = 1 only the counter-rotating 
half and for i^ = the full possible range of angular momenta is populated. Components with = are therefore 
non-rotating. An impression of the appearance of those families of components can be obtained by considering their 
parameters one by one. The parameter a indicates the degree of central concentration and /3 the degree of rotation in 
both families. For non zero /3, the density distributions become toroidal. For the second family (F2) zo is the absolute 
vertical cut-off for the component and 7 determines the vertical scaleheight. The larger 7, the faster the density falls 
off with increasing height above the plane. A smooth transition for p — >0 at zq for all R is ensured by the functional 
form of the second family. 

The program QP reads the allowed values for all parameters from an input library. For each component it determines 
the coefficient C for which the value of D (Eq. 2) is minimized. The component with smallest D is then chosen as the 
first in the series that forms the total distribution function. Subsequently, all the remaining components are checked 
for the smallest value of D in combination with the first component. The coefficient of the first component does 
not have to remain fixed; it can even become zero in the process of converging. The sum of the first N components 
has to be positive for all TV, so that the series gives a valid distribution function - positive everywhere in phase 
space - at any instant in the convergence, but individual C's can be negative in principle. In fact, one may use the 
negative coefficients to test if the model components have any physical meaning; when positive coefiicients alternate 
with negative ones in successive components, the model is similar to a power-series development and thus a purely 
mathematic construct. This means that none of the components indivually match the data well and one may want to 
consider a new component library. We demand that all coefiicients be positive for all our models, but test the final 
library allowing also negative coefficients. 

The best solution, for given input library, is reached when the value of D (Eq. 2) has converged to within a few 
percent. The number of components in the converged solution is mostly of the order of = 5 for the models in this 
paper. We ran QP with a great variety of input libraries, starting with one that spans a wide range for the parameters 
and fine tuning toward preferred solutions. For example, if a = 20 is selected from a library with a = (3, 10, 20, 30) 
then the next library will have a = (15,20,25). These libraries are relatively small, with ~ 60 components, to 
reduce the computing time (t cx Ncl/{Nc — Mc)!). One should be careful not to exclude non-preferred values for the 
parameters once and for all, because QP seeks out the best combination of components. It may be that in the first 
trial runs, for example, F2 components with a > 10 are never used. After the fine-tuning process, such components 
could nevertheless improve the solution when combined with the final components. We therefore always ran QP with 
a large library {Nc ~ 250), that combined the best library with earlier ones, as a final test. Mostly, the solution in 
these final runs did not differ from that obtained with the best (small) library. 

In Fig. 4 we show how orbits populate different regions of phase space. The energies and angiilar momenta that 
an orbit can have are determined by the potential. A detailed explanation is given in the figure caption. 
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Fig. 4. Orbits are points in the (i?, L^) or phase space. For the AX2 potential they are arranged as follows. By 
convention, co-rotating orbits have negative angular momentum. The outer thin line indicates the {E, L^) for circular 
orbits; the shape is determined by the potential and thus only valid for AX2. Only the regions inside this line can be 
populated, as circular orbits have maximum Lz for given E. The dashed ellipses are the perimeters of regions accessible 
for orbits that pass through the plane at a given radius R\ their intersections with the circular-orbit line obviously 
indicates the E, L.^ of the circular orbit at R. The dot dashed ellipses are the perimeters of regions accessible for orbits 
that pass through the z = A kpc plane at given radius R (same i?'s as for z = kpc). The thin line reminiscent of 
the circular-orbit line connects the maximum for orbits with different E, that reach z = 4 kpc somewhere along 
its orbit. For =0, the orbits are radial (vertical dashed line). Isotropic distributions {(3 = 0) are independent of L^, 
and will give horizontal contours in this diagram. Rotating, disky distributions yield contours roughly parallel to co- 
and/or counter rotating half of the circular-orbit line. If the distribution is flattened, say z < 4 kpc, the phase-space 
density will only be non zero between the two thin, solid lines (for maximum angular-momentum at z = and 2 = 4, 
respectively). is scaled to increase the relative resolution in the inner regions (£^ ~ 1, Lz ~ 0). 

The components arc integrated out to a predefined limit. For the models using the whole AOSP sample we use a 
horizon at 13 kpc; in every direction the model is integrated out to 13 kpc from the position of the observer (Fig. 5). 
This limit is chosen somewhat larger than the observational limit (Sevenster 1999). All galactic radii smaller than 5 
kpc are thus sampled twice at each line of sight at \£\ < 39° and best constrained; radii between 6 kpc and 8 kpc are 
sampled at once or twice per line of sight (Fig. 5). Radii larger than 8 kpc are sampled only once for lines of sight at 
1^1 > 36°. At the largest longitude used in the modelling, 45°, the horizon lies at a galactic radius 9.2 kpc. In models 
for only low-outflow sources we use a horizon at 11 kpc, as these have a smaller observation limit (Sevenster 1999). 

4. Results 

The best 21 distribution function (DFA) was obtained using the component library given in Table 4. Table 5 gives 
the components of DFA with their coefHcients, C, as well as the masses, M^, of the components in the region R <8 
kpc and \z\ <4 kpc . The total DFA is the sum of Cn-Fn, but the actual relative contribution of component to 
DFA is Cn-^w,n- In Fig- 6 we show the phase-space density of DFA. The combined projected moments (Sect. 3) of 




Fig. 5. The hatched wedge gives the part of the plane over which the distribution function is integrated to yield the 
moments, mimicing the observations with a horizon at 13 kpc (dashed circle d=13 kpc, see Sect. 3.3). The galactic 
Centre is at the intersection of the straight lines. Evidently, not all galactic radii in the observed region are sampled 
equally well. Radii < 5 kpc are sampled optimally, with two intersections at any longitude (or four for < 1.4 kpc). 
The largest radius that is sampled is 9.2 kpc. 



Table 4. The input library. Explanation see Sect. 3.3. 



Family 


a 


/3 


7 


20 


is 


iVc 


Fl 


5, 10,20 


0,1,2 






0,-1,1 


27 


F2 


2,4, 8, 10 


1,2 


2,6 


1,5 


-1 


32 



Table 5. The components (in order of choice by QP from Table 4) of the best-fit distribution function DFA (Sect. 4) 
for a horizon at 13 kpc. C gives the coefficient of the component, its mass within a fixed cylinder (see Sect. 4). 



Family 


a 


/3 


7 


20 


is 


c 


Mw 


F2 


2 


2 


6 


1 


-1 


1.94E4 


3.35E-2 


Fl 


10 












3.52E2 


5.69E-1 


Fl 


20 


1 






-1 


1.54E6 


2.89E-5 


Fl 


20 












5.92E3 


9.69E-3 


Fl 


5 












1.43E1 


1.24E+1 



DFA and of the data are shown in Fig. 7 and the true projected moments (as used in the fit) in Fig. Al (Appendix 
A). At the inclusion of the fifth component, the value of D (Eq. 2) has converged to within 2%, to 22% of the initial 
value. As explained in Sect. 3.3, we did not allow negative coefficients for the components in the model distribution 
function. We tested that the outcome (Table 5) is not dependent upon this; exactly the same results are obtained 
when negative coefficients are allowed. The small-scale (non-axisymmetric) features are, correctly, mostly neglected 
by QP. Apart from this, the main discrepancies between data and model are seen in the scaleheight (Fig. 7c, k), the 
central dispersion (Fig. 7g,h) and the vertical rotation profile at \i\ ~ 8° (Fig. 7f). 

The underlying reason is the same for all these discrepancies: the dispersion is too high to be explained by a 21 
model that fits the other moments. The line of sight to the galactic Centre is parallel to the radial direction, so the 
observed central dispersion ctq depends on the radial dispersion ctr only. Since the scaleheight ft-z oc a\ and in 21 
distributions ctr = ctz, a component that increases ctq will increase inevitably, via ctr and Uz, the scaleheight at ^ = 
0°. The bad reproduction of the vertical kinematic profiles arises because components that give cylindrical rotation 
(mainly F2) have low dispersion. Components that give vertically-constant dispersion profiles (mainly Fl), as observed 
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L,«V(2E/E„) 

Fig. 6. The phase space density as function of E, L.^ (see Fig. 4) for DFA. the best fit two integral model. The 
logarithmic contours range from 1 (top) to 10^ (bottom, arbitrary units), so the highest phase-space density is at 
{E, iz;)=(0. 0,1.0) : the central region of the Galaxy. Comparing to Fig. 4; one can clearly discern a thin disk (closely- 
spaced contours at negative angular momentum) and an isotropic component (horizontal contours). 



for 1^1 > 15°, do not have cylindrical rotation. Note that the flatness, high central dispersion and cylindrical rotation 
are all signs of the barred central Galaxy (eg. Kormcndy 1993). 

For oversmoothed data (with a kernel twice as large as the optimal kernel discussed in Sect. 3.1), the surface- 
density vertical profiles are fitted much better (Fig. A2c,k,l). Accordingly, the model cto is indeed higher (Fig. A2g,h) 
and the model rotation is more cylindrical (Fig. A2f). The fit to the dispersion is better as well, because the central 
dispersion is no longer so sharply peaked. 

In Fig. A3 we present the model derived for the data smoothed with elongated kernels (see Sect. 3.1). The central 
scalehcight is now very small and the vertical surface-density profile is only fitted above I6| = 1°. Therefore, the 
central dispersion could be fitted well; it is also lower than in the standard-smoothed data because the (central) 
disk contributes more, decreasing the dispersion in the plane. For |6| > 1°, the minor-axis surface-density profile is 
modelled well at the cost of the modelled minor-axis dispersion. In Fig. A4 we show the best model using the BD2 
potential (Sect. 3.2). The global rotation and dispersion are not fitted well, as expected from the discrepancy between 
potential and data (Fig. 2). 

4.1. Errors, biases and stability 

To estimate the errors in the input data and in the distribution function, we applied the "bootstrap" method (Press 
et al. 1992). New samples were created by drawing randomly 507 stars from the data sample (consisting of 507 stars) 
"with replacement" . This means that the same star can appear in the sample more than once. We are allowed to do this, 
because the sample is virtually free of any biases (S97A, S97B). These 25 new samples were smoothed and modelled 
in exactly the same way as the original data. The mean of and the scatter in the results provide biases and error 
bars on the input gridded data as well as on the model distribution function and its moments (Fig. 7,Al): bootstrap 
measurements have the same distribution with respect to the original measurement as the original measurement has 
with respect to the "true" value. The errors arising from the gridding are small and, as expected, the data are free of 
bias, except for some small latitude dependent bias (S97A, S97B; most notable in Fig. 7f). For the model, the small 
errors and biases indicate that the distribution function is in general well constrained by the data and close to the 
"true" distribution function, given the limitation of being a two-integral model. 

Some of the model moments are considerably biased, though (Fig. 7, Al). The mean vertical profiles in Fig. 7j,l 
are actually closer to the data than the best-fit model. This means that these deviations between the DFA moments 
and the data are not due to an intrinsic limitation of our modelling technique, but that for some reason the data cause 
the model to be biased, eg. via the inevitable limitations of observational sampling. One should be careful not to think 
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Fig. 7. Cuts through combined projected moments of the model distribution function DFA (dashed) and through the 
data (sohd; see Sect. 4). The means and error bars on data (triangles) and model (circles) are found via bootstrapping 
(note that the curves are not the means, see Sect. 4.1). The data was adaptive-kernel smoothed with initial kernels of 
l°xl°x 30kms~^ (see Sect. 3.1). Labels in the panels indicate for what value of latitude or longitude, respectively. 
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The (a;, y) and {x, z) positions for 1000 particles from the random realization of DFA. 
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Fig. 9. The cumulative-number density as a function of R and of z for the initial distribution of the random realization 
(solid curve) and the distribution after ten galactic revolutions (dashed curve). According to the Kolmogorov-Smirnov 
test, the probability that initial and final distributions are the same is more than 95%. 



that the "average model" is closer to the "true" distribution function than DFA is; it only means that DFA is not 
close to the "true" distribution function for \£\ J; 35°. 

On the other hand, the biases on the model in Fig. 7c, i are away from the data. Hence the modelling technique 
really cannot provide good results here; one might guess using only two integrals is the inhibiting factor here. There 
is also a small bias (^ 1 a) away from the data on the longitude model profiles in Fig. 7a, b ; this could be due to the 
obvious non-axisymmetric features present in the data and thus the limitation of our axisymmetric model. Similarly, 
the bias at larger longitudes on the model dispersion profiles (Fig. 7g,h) could be caused by the limitations of our 
model potential that is not quite adequate at larger radii. So, although globally the 21 model reproduces the observed 
moments well and with small uncertainties, it is presented with real problems by the vertical profiles at £=0° (Fig. 7c, i), 
as well as by the cylindrical rotation at \(.\ = 8° (Fig. 7f). 

We created an N-particlc realization (N=5GGG) of the distribution function (see eg. van der Marel et al. 1997) and 
evolved it, in the global potential used in QP, to show that DFA is numerically stable, as it should be as a valid 
function of the integrals of motion. The realization is shown in Fig. 8, where we plot 1000 particles in the {x, y) plane 
and in the (x, z) plane. As criteria for stability we checked the total radial and vertical cumulative-density profiles and 
the total energy. To quantify the stability of the density profiles we use the Kolmogorov-Smirnov (KS) test described 
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Table 6. The components (in order of choice by QP) of the best-fit distribution function for the galactic-centre 
sample, integrated over a spheroid of 0.5x0.25 kpc centered at the GC. 
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by Press et al. (1992). The radial and vertical profiles pass the KS test very well: the differences between the initial 
distribution and that after ten galactic revolutions are entirely negligible (Fig. 9, KS probability ^95%). The total 
energy shows no variations other than of the order of the accuracy of the integration (10~^). 

5. Third integral 

The problems of the 21 fit, mentioned in Sect. 4, may be overcome by the use of three-integral (31) models, as for these 
the radial and vertical dispersions do not have to be the same and therefore the radial and vertical distributions are not 
coupled. For axisymmetric Stackel potentials three integrals of motion - E, I2 (= O.SX^) a-nd I3 - are known analytically 
(de Zeeuw 1985). We use QP with three types of 31 orbital components; two are derived from the components in Eq. 3,4. 
The third-type components are axisymmetric Abel components (Dejonghe & Laurent 1991). Using the same input 
data as for DFA we obtain the distribution shown in Fig. 10. The 31 distribution function comprises a few components 
with small, negative coefficients. The problems of the fit in DFA (Fig. 7) are largely solved; in the inner regions, 
the dispersion is 20 kms~^ higher for all latitudes and at the same time the minor-axis surface-density profile is 
fitted better. Also the cylindrical rotation at intermediate longitudes (Fig. lOf vs. Fig. 7f) is reproduced better; the 
deviation of the model rotation from the data at 6 = 3° is only half that in the 21 model. The components of the 31 
distribution function are only truly 31 inside ~ 4 kpc; outside that the dependency on the third integral decreases 
until it disappears at the solar radius. It is well known, however, that also in the Disk a third integral is needed to 
describe the distribution of most populations, because the local vertical dispersion does not equal the radial dispersion 
(eg. Wielen 1977). We will not draw conclusions about the third integral in the Disk, as the AX2 potential does not 
represent properly the gravitational forces outside radii of ^5 kpc. Clearly, in the innc;r regions a third integral is 
needed to model the observed high dispersion and small scaleheight. An upcoming counterpart to the AOSP survey 
covering positive longitudes (see S97A,B) will enable us to construct a proper triaxial, three-integral model. The 
observed positive negative longitude asymmetries in stellar kinematics, in combination with the known asymmetries 
in the stellar surface density, will be essential to do this. 

6. The central 100 pc 

In Sevenster et al. (1995; SDH) a sample of 134 stars in approximately the inner square degree of the Galaxy (Lindqvist 

et al. 1992) was modelled using the BD2 potential (Table 2), with an additional Plummer potential truncated at 100 
pc. With the AX2 potential we obtain virtually the same results, this time without an additional Plummer potential. 
For the same data gridding as in SDH, over 20 nearest neighbours, there are co- and counter-rotating F2 components 
of similar extent (ie. similar a) and more concentrated co-rotating or fully isotropic Fl components. 

It is difficult to determine the correct way to smooth this small sample and the exact results are dependent upon 
the smoothing. The model we give in Fig. 11 (Table 6) is obtained from the data averaged over 30 nearest neighbours. 
In this case no counter- rotating components are found by QP. The value of D (Eq. 2) has converged to within 1% at 
the inclusion of the fifth component (Table 6), but the final value of D is 60% of the initial. The latter indicates that 
this fit can be improved upon; there are only components with a = 200 and a = 900 in the input library so this may 
have been crude. The radial scale of the a = 200 components is roughly 200 pc; that of the a = 900 components 50 
pc. The rotation in the inner 100 pc comes from the Fl components. 

This is even more pronounced if we use only the 10 nearest neighbours to calculate the moments. In this case the 
observed dispersion at b=0° increases from 50 kms"'^ at £ = 0° to 120 kms~^ at £ = 0°.6, which is the observed 
value for the AOSP sample (Fig. 7). The a = 200 are co- and counter-rotating, respectively, as in found in other runs, 
mainly to fit this increasing dispersion. All rotation comes from a = 900 Fl components in this case. 

Possibly the a = 200 components are connected to a part of the GC sample that forms the innermost extension 
of the Bulge; the kinematics, radial scale and large vertical scale {zq = 1 kpc) fit in well with that. The counter- 
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Fig. 11. The fit to the sample of galactic centre stars (data averaged over 30 nearest neighbours; see Sect. 6). The 
cuts are at different longitudes and latitudes than in Fig. 7. 



rotation in these components found for a range of inputs is used to increase the dispersion artificially; possibly the 
central concentration of the potential is not fully adequate; or a third integral, that reproduced the large AOSP central 
dispersion, is needed also in this model. The a = 900 components may form the true GC sample (SDH; Sjouwerman et 
al. 1998a, b; Sevenster 1999), judging from their high rotation and small radial and vertical scales. These components 
are probably formed mainly by high-outflow sources. 



7. Discussion 

7.1. Density 

For the 21 as well as the 31 model, the scalelength and scaleheight at i? = 6 kpc are 2.5 kpc and 200-250 pc, respectively. 
At i? = kpc, they are 200 220 pc and 150 pc. These scales arc very similar to those foimd for the same sample in 
an analysis of the surface density only (Sevenster 1999), except for the scalelength that is smaller but still large with 
respect to most determinations of the scalelength (see Sackett 1997), although Binney, Gerhard & Spergel (1997) also 
find a scalelength of 2.5 kpc from fitting COBE data. 

7.2. Orbits 

The fractions of co-rotating and of eccentric orbits, respectively, can be defined by the following formulae : 



f{E,L,)dL, / / f{E,L,)dL, 

"i'max / *' "-^Tnax 

f{E,L,)dL,/ / f{E,L,)dL, 



F = 

^ ecc — 
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Fig. 12. (a) The fractions of the phase space density on co rotating orbits (Eq. 5, thick sohd curve) and on eccentric 
orbits (Eq. 6, thick dashed curve). The sum of the two (thin sohd curve) equals 1 only if the distribution is fully isotropic, 
fully co-rotating or contains only eccentric orbits. The biases (points) and error bars are obtained via bootstrapping 
(Fig. 7, Sect. 4); the thin dashed curve gives the sum of the two bootstrap means. DFA is fully isotropic in the centre 
{E = 1), but the bias indicates that the AOSP sample is not. (b) The mean rotation for DFA and the 31 model, as 
well as for the Fl and F2 components of DFA separately, is shown as a function of radius. 



where Lmax is the absolute value of the angular momentum of a circular orbit with energy E. In Fig. 12a these 
fractions are shown for DFA. In the centre {E/Eq > 0.95, which coincides for circular orbits with i?cir < 160 pc, for 
radial orbits with i?rad < 230 pc. Fig. 4) the distribution function is isotropic to within 1%. For E/Eq < 0.5 ( i?cir > 
2.5 kpc, i?rad > 4 kpc) more than 99% of the mass is on almost-circular, co-rotating orbits. The biases and errorbars 
indicate that the distribution function at high binding energies is not very well-constrained. This can be seen in the 
moments only for the surface density, for which the errorbars near the galactic Centre are large (Fig. 7a,b,c). 

There is considerable bias on the fraction of DFA on eccentric orbits, which means that this fraction is not close 
to the "true" fraction of eccentric orbits in the central Galaxy. The fraction of co-rotating orbits is less biased and 
always larger than 50% (in the mean), so there is no significant net counter rotation in the AOSP sample. 

There is a turn-over in the energy distributions in Fig. 12a at E/Eq ~ 0.7 (i?cir= 840 pc, i?rad= 1-4 kpc) and 
the fractions of eccentric orbits and of counter rotating orbits decrease quickly outside this radius. Around R = 800 
pc, the rotation curve of the Fl component (Fig. 12b) reaches its maximum (31.5 kms~^ ). The rotation is continued 
fairly smoothly by the F2 component, though. A disky and a bulge like regime can be identified in energy, but not so 
clearly in radius, although a mild transition can be seen at a radius 2.5 kpc (Fig. 9). There is an "isotropic-rotator" 
regime inside ~ 1 kpc outside which the disk starts. Between 1 kpc and 4 kpc the isotropic components contribute 
to a non-negligible fraction of mass on eccentric and counter-rotating orbits. The rotation is almost linear out to 2.5 
kpc, outside 4 kpc the regime is purely disky. This is in agreement with arguments that barred bulges do not extend 
beyond their co-rotation radius (eg. Contopoulos & Grosb0l 1986; Elmegreen & Elmegreen 1985), which in the Galaxy 
is at 4 to 5 kpc. 



7.3. Dispersions 

The dispersions that derive from DFA are obviously biased by the fact that o"r=(Tz. At the solar radius, DFA yields 
(cR,cr0,crz) = (11 kms^^, 22 kms^^ , 11 kms^-^) so (7p = y^cr^ + '^r) = 25 km s^-^ (Fig. 13a). This cTp agrees 
with a population of ~ 1.5 Gyr for which the observed full velocity ellipsoid is (19 kms~^ , 15 kms~^ , 10 kms~^ ) 
(Wielen 1977). So is forced to reproduce most of trios as ctr cannot be larger than a^, that in turn is limited by the 
scalcheight. In the 31 model, the local dispersions are all 13 kms"-*^ . For yoimg stars (< 0.5 Gyr), the three dispersions 
are observed to be almost equal, but have a much lower value. The 31 dispersions yield CTp = 18 kms~^ , equivalent 
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Fig. 13. The radial (dashed), the azimuthal (solid) and the vertical (dotted) dispersions in kms ^ at 2; = for DFA 
(a) and the 31 model (b) as a function of radius. For DFA the vertical dispersion equals the radial. 



Table 7. The components (in order of choice by QP from Table 4) of the best-fit model to the high-outflow sources 
(horizon 13 kpc). 
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to observations of stars of ~ 1 Gyr, for which the full velocity ellipsoid is (14 kms~^ , 11 kms~^ , 8 kms~^ ) (Wielen 
1977). Despite the decoupling from a^, ctr is still not larger than a^. As we already noted in Sect. 5, the 31 model does 
not contain 31 components at radii larger than ^ 4 kpc. Although individually ctr and cr^ do not match observations, 
the dispersion in the plane CTp has the value expected for a population of the average age of the OH/IR stars in the 
Disk (^ 1.5 Gyr, Sevenster 1999). Keeping in mind that our model potential is not fully adequate at radii larger than 
5 kpc, the "deviant" model-dispersion ratios may indicate that a significant number of OH/IR stars is not on epicyclic 
orbits. 

A well-observed dispersion is the line-of-sight dispersion toward Baade's window {£ = l°,b = —4°) of 113tg 
kms"""^ (Sharpies ct al. 1990). DFA yields 107 kms"-'^ and the 31 model 113 kms"-*^ , hence the observed dispersion, 
even at the higher latitude of Baade's window, is matched somewhat better by the 31 model. The proper-motion 
dispersions toward Baade's window are (cr^, ab) = (3.2±0.1 masyr~^,2.8±0.1 masyr"^) (Spaenhauer et al. 1992). DFA 
and the 31 model yield (1.6,1.5) and (3.5,2.4), respectively (taking the detectability of sources inversely proportional 
to the distance squared and integrating out to 8 kpc). 

7.4- Disk versus Bulge 

Imagine that the total distribution function, DFA, indeed fully describes the stellar dynamics of the Galaxy. Each 
OH/IR star can be thought of as a random realization of DFA - more specifically, as drawn from one of the DFA 
components. The probability, then, that a star S is drawn from, say, component DFAl, is the conditional probability 
(SgDFA1|SgDFA). This, according to Bayes' rule, is proportional to the density of the component at the position of 
S. Hence, the component with the highest density, integrated over the unmeasured coordinates, at the position of a 
star, is most likely to have "generated" that star. 

In Fig. 14 we show this by plotting for each component of Table 5 the longitude-latitude diagram and longitude- 
velocity diagram of the stars for which this component gives greatest probability. Although one should be careful to 
interpret the components of the distribution function exactly as physical components of the Galaxy, it is clear that 
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Table 8. The components (in order of choice by QP from Table 4) of the best-fit model to the low-outflow sources 
(horizon 11 kpc). 
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the first component DFAl forms the main galactic Disk and DFA2-4 are connected with a slowly rotating isotropic 
bulge. The role of DFA5 is difficult to assess. The fraction of stars with longitudes below —20° that connect to DFA5 
is 25%. However, the mass fraction of this component is of the order of 2% for radii larger than 2.5 kpc (Fig. 12). The 
Disk stars connected to DFA5 are probably mainly those that make up the local features that are not fitted by DFA 
(Fig. Al). As those stars have rather deviant velocities, they fit in best with DFA5, because it has very high velocity 
dispersion (~100 kms^^ ). This does not mean that those stars instigated the inclusion of DFA5 in DFA, after all 
they are not properly represented by the flt. DFA5 may represent the tail of the Bulge, that apparently protrudes to 
R ^4: kpc (cf. Fig. 12a). Indeed the stars connected to DFA5 have the lowest total probability to be connected to 
DFA, supporting the idea that they were not really fitted very well. The most probable component is DFA4, followed 
by DFAl. The star with the lowest probability to come from DFA is the one in the extreme lower right corner of 
Fig. 14e,f,k,l . 

We tentatively connect the Fl components to the Bulge and the F2 component to the Disk. As DFAl is well 
constrained by stars in regions where only the Disk is contributing, we may conclude that also at lower longitudes it 
represents the (foreground) Disk, as is supported by Fig. 14a, g. This means that 50% of the AOSP sample is identifled 
with the Disk and that the fraction of Disk stars at \£\ <5° is 20%. 

7.5. Young versus Old 

In Fig. 14 different symbols are used for high-outflow, low-outflow and single-peaked sources, respectively. OH/IR 

stars with high outflow velocities are in general younger than those with low outflow velocities (see Sevenster 1999 
and references therein). The single-peaked sources can be either (very) young or (very) old; we will not include them 
in this discussion. All Fl components of DFA have more low-outflow- than high-outflow sources (4:3) connected to 
them. For the F2 component this is exactly the other way around, as expected from the fact that it connects to Disk 
sources mainly. We can assess further the age-dependence of the distribution function by modelling the two groups, 
separated in outflow velocity at 14 kms~^ (excluding single-peaked sources), individually, with the same potential 
and library we used to obtain DFA. The results are given in Table 7 and Table 8, respectively. 

The main disk (DFAl) and the rotating-bulge component DFA3 return for both the low-outflow- and the high- 
outflow sources. DFA4 returns for the high-outflow sources only and DFA2 for the low-outflow sources. These com- 
ponents are the main Bulge components, a younger and an older (more extended), respectively, forming the "isotropic 
rotator" together with DFA3 (Sect. 7.2). The high-outflow sources have an extra F2 component with the same vertical 
extent as the main DFA disk, but more centrally concentrated and less strongly rotating. This may be the youngest 
part of the Bulge that is not massive enough to be seen in DFA. The low-outflow sources have an extra F2 component 
that has the same radial extent and rotation (/?) as the main DFA disk, but is less concentrated toward the plane 
(/iz=300 pc). This may be the older Disk, heated from the flatter Disk. Again it is not massive enough to appear in 
DFA. 

In Fig. A5&A6 we show the fits of DFA and the distribution function of Table 8, respectively, to a sample of OH/IR 
stars that reaches higher latitudes than the AOSP sample but is incomplete in the plane (|6| < 3°; te Lintel Hekkert 
et al. 1991). The coefficients for the components of DFA and Table 8 are redetermined for this sample. We use the 
horizon that optimizes the fit to the surface density (13 kpc). Clearly, a thicker disk component such as seen in the 
low-outflow sources is essential to explain the still cylindrical rotation in the Lintel sample; in fact the coefficient for 
the F2 component with 7 = 6 is zero for the fit in Fig. A6 (as this flat component would be severely undersampled 
by the Lintel sample). At the higher latitudes of the Lintel sample the stars are on average older, like the low-outflow 
AOSP sources, and these older stars apparently need a thicker disk component, primarily to describe their kinematics. 
The transition between the Disk and the Bulge (Fig. 12), at least at higher latitudes, is more continuous than suggested 
by DFA (Fig. A5e vs. Fig. A6e; see discussion in Sect. 7.2). 
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Fig. 14. This figure shows, for each component (from the top down in the order of Table 5; in the bottom panel for 
DFA in total), the longitude-latitude diagram (left) and longitude-velocity diagram (right) of the stars "associated" 
with it. The different symbols indicate high-outflow (triangles), low-outflow (circles) and single-peaked (crosses) 
OH/IR stars. For further explanation see Sect. 7.4. 



In summary, it seems that there are no distinct dynamical components such as a Bulge or a thick Disk. In fact, 
the young Bulge has a part that in vertical extent is very similar to the Disk, in radial extent to the more isotropic 
older Bulge and intermediate in its degree of rotation. The older Disk is very similar in radial extent and rotation to 
the younger disk, only a little thicker and as such connecting even more smoothly to the Biilge. We do not sample the 
very old Bulge (^10 Gyr), that may be the inner halo or 'V^-*-^ spheroid" and was found to be dynamically different 
from the younger "nuclear Bulge" by Rich(1990). 

The connection we find between Bulge and Disk, especially their similar vertical extent for the younger stars, is 
in agreement with the notion that the Bulge is triaxial and that this Bar formed via disk instability (see Sevenster 
1999). Of course, we have already seen several signatures of the existence of the Bar, in the need for a third integral 
to explain the dynamics in the central degrees and the the cylindrical rotation (sec discussion in Kormcndy 1993). 

Note that in this case, the fraction of foreground disk stars we estimated earlier for the central 10° (20%) may be 
too high, as part of the disk-like Bar stars would probably be modelled by DFAl as Disk stars. 
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8. Conclusions 

Using a simple, axisymmctric potential we construct a stable two integral distribution function (DFA) that gives a 
very good global fit, in the first three projected moments, to our "AOSP" sample of OH/IR stars in the plane. Some 
detailed discrepancies between the model DFA and the data indicate that the distribution of OH/IR stars is influenced 
by the barred potential in the inner regions of the Galaxy. A three-integral model improves the fit for the inner regions 
considerably, even with the same axisymmetric potential. Durand et al. (1996) also concluded there is a need for a 
third integral, from similar work on distribution functions, using planetary nebulae. 

The energies of stars in the plane seem to separate into a bulge-like and a disky regime at E/Eq ^ 0.7. This 
separation, as seen in DFA, is too distinct, however, when compared to the kinematics of an OH/IR star sample at 
higher latitudes. We conclude there is no evidence for discrete large components in the inner plane. On the contrary, 
models of several subsamplcs of younger and older OH/IR stars suggest that the Disk and the Bulge are very similar. 
We confirm the result of Sevenster et al. (1995) that a sample of galactic-centre OH/IR stars may consist of the 
inner-most part of the Bulge plus an extra component. The latter is the only truly distinct dynamical component in 
the inner galactic plane. 
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Appendix A Figures for derived models 

In this appendix we show the figures of the cuts in longitude and latitude for a variety of models discussed in the 
main text. 
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Fig.A2As Fig. 7, for best-fit distribution function for data smoothed with spatial initial kernel of 2° (same library as DFA; 
Table 4). 
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Fig. A3 As Fig. 7, for best-fit distribution function for data smoothed with spatial initial kernel of 1°:0.5° (same library as 
DFA; Table 4). 



22 



Sevenster et al.: Distribution functions for evolved stars 




"> - e = 8.00 



1.5 2 2.5 3 (k) 0.5 

:,jd. (■) 



1.5 2 2.5 3 (I) 0.5 

ud. (■) 



bdok7_0909_215a 
1.5 2 2.5 3 



Fig.A4As Fig. 7, for best-fit distribution function with BD2 potential (Table 2) and the combined component library (258 
components, Sect. 3.3). 
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Fig. A5 Fit of DFA (Table 5, different coefficients) to the Lintel sample. The horizon is at 13 kpc. 
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Fig. A6 Fit of the low-outflow-AOSP distribution function (Table 8, different coefficients) to the Lintel sample. The horizon is 
at 13 kpc. 



